home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-04 / amsf20.zip / AMST3.FOR < prev    next >
Text File  |  1992-01-06  |  7KB  |  233 lines
  1. C     ******************************************************************
  2. C     *                                                                *
  3. C     *                         E I G E N                              *
  4. C     *                                                                *
  5. C     ******************************************************************
  6.       PROGRAM EIGEN
  7.       IMPLICIT INTEGER*4(I-N)
  8.       IMPLICIT REAL*8 (A-H,O-Z)
  9.       COMMON  MAVAIL,IA(30000)
  10.       MAVAIL=30000
  11. C ... SAMPLE: EIGEN PROBLEM SOLVER
  12.       PRINT *,' '
  13.       PRINT *,'EIGEN PROBLEM SOLVER'
  14.       PRINT *,' '
  15.       PRINT *,'ENTER EQUATION NUMBER '
  16.       READ  *, N
  17.       PRINT *,'NOW ENTER STIFFNESS MATRIX & MASS MATRIX'
  18.       PRINT *,' '
  19. C ... TEST IN-CORE DATA MANAGEMENT
  20.       CALL DBOPEN(1,'TSTDAT.DAT','NEW')
  21.       PRINT *,'DATABASE OPENED.'
  22.       CALL DEFINE(1,'STIF',0,1,N,N,0,NA)
  23.       CALL DEFINE(1,'MASS',0,1,N,N,0,NB)
  24.       CALL DEFINE(1,'VCTR',0,1,N,N,0,NC)
  25.       CALL DEFINE(1,'EIGN',0,1,N,1,0,ND)
  26.       CALL DEFINE(1,'WORK',0,1,N,1,0,NE)
  27.       print *,'DEFINE MATRICES DONE.'
  28.       CALL MATINP(1,'STIF')
  29.       CALL MATINP(1,'MASS')
  30.       CALL MATOUT(1,'STIF')
  31.       CALL MATOUT(1,'MASS')
  32.       RTOL = 1.0D-12
  33.       CALL JACOBI(IA(NA),IA(NB),IA(NC),IA(ND),IA(NE),N,RTOL,15,0,6)
  34.       CALL MATOUT(1,'EIGN')
  35.       CALL MATOUT(1,'VCTR')
  36.       CALL DIR(0)
  37.       CALL DBCLOS(1,'DELETE')
  38.       STOP 'DONE.'
  39.       END
  40. C     ******************************************************************
  41. C     *                                                                *
  42. C     *                          J A C O B I                           *
  43. C     *                                                                *
  44. C     ******************************************************************
  45.       SUBROUTINE JACOBI (A,B,X,EIGV,D,N,RTOL,NSMAX,IFPR,IOUT)
  46.       IMPLICIT INTEGER*4(I-N)
  47.       IMPLICIT REAL*8 (A-H,O-Z)
  48. C
  49.       DIMENSION A(N,N),B(N,N),X(N,N),EIGV(N),D(N)
  50. C
  51. C     TO SOLVE THE GENERALIZED EIGENPROBLEM USING THE GENERALIZED
  52. C     JACOBI ITERATION
  53. C
  54. C     INPUT VARIABLES :
  55. C       A(N,N)   = STIFFNESS MATRIX (ASSUME POSITIVE DEFINITE)
  56. C       B(N,N)   = MASS MATRIX      (ASSUME POSITIVE DEFINITE)
  57. C       D(N)     = WORKING VAECTOR
  58. C       N        = ORDER OF MATRIX A AND B
  59. C       RTOL     = CONVERGENCE TOLERANCE (SET TO 1.E-12)
  60. C       NSMAX    = MAXIMUM NUMBER OF SWEEP ALLOWED (SET TO 15)
  61. C       IFPR     = FLAG FOR PRINTING DURING ITERATION
  62. C                  ( EQ.0 NO PRINTING, EQ.1 PRINT INTERMEDIATE RESULTS)
  63. C       IOUT     = OUTPUT LUN
  64. C     OUTPUT VARIABLES :
  65. C       A(N,N)   = DIAGONALIZED STIFFNESS MATRIX
  66. C       B(N,N)   = DIAGONALIZED MASS      MATRIX
  67. C       X(N,N)   = EIGENVECTORS STORED COLUMNWISE
  68. C       EIGV(N)  = EIGENVALUES
  69. C
  70.  
  71. C     INITIALIZE EIGENVALUE AND EIGENVECTOR MATRICES
  72.       DO 10 I=1,N
  73.       IF (A(I,I).LE.0..OR.B(I,I).LT.0.) THEN
  74.          WRITE(IOUT,2020)
  75.          STOP 'JACOBI'
  76.       END IF
  77.       D(I)=A(I,I)/B(I,I)
  78.  10   EIGV(I)=D(I)
  79.       DO 30 I=1,N
  80.       DO 20 J=1,N
  81.  20   X(I,J)=0.
  82.  30   X(I,I)=1.0
  83.       IF(N.EQ.1) RETURN
  84. C     INITIALIZE SWEEP COUNTER AND EIGEN ITERATION
  85.       NSWEEP=0
  86.       NR=N-1
  87. C
  88. C     WE START ITERATION
  89.  40   NSWEEP=NSWEEP+1
  90.       IF(IFPR.EQ.1) WRITE(IOUT,1000) NSWEEP
  91. C     CHECK IF PRESENT OFF-DIAGONAL ELEMENT IS LARGE ENOUGH TO
  92. C     REQUIRE ZEROING
  93.       EPS=(0.01**NSWEEP)**2
  94.       DO 210 J=1,NR
  95.       JJ=J+1
  96.       DO 210 K=JJ,N
  97.       TT=A(J,K)*A(J,K)
  98.       TB=A(J,J)*A(K,K)
  99.       EPTOLA=DABS(TT/TB)
  100.       TT=B(J,K)*B(J,K)
  101.       TB=B(J,J)*B(K,K)
  102.       EPTOLB=TT/TB
  103.       IF ((EPTOLA.LT.EPS).AND.(EPTOLB.LT.EPS)) GO TO 210
  104. C     IF ZEROING IS REQUIRED, CALCULATE THE ROTATION MATRIX ELEMENT CA, CG
  105.       AKK=A(K,K)*B(J,K)-B(K,K)*A(J,K)
  106.       AJJ=A(J,J)*B(J,K)-B(J,J)*A(J,K)
  107.       AB=A(J,J)*B(K,K)-A(K,K)*B(J,J)
  108.       CHECK=(AB*AB+4.0*AKK*AJJ)/4.0
  109.       IF (CHECK) 60,60,70
  110.  60   WRITE (IOUT,1004) CHECK
  111.       STOP
  112.  70   SQCH=DSQRT(CHECK)
  113.       D1=AB/2.0+SQCH
  114.       D2=AB/2.0-SQCH
  115.       DEN=D1
  116.       IF (DABS(D2).GT.DABS(D1)) DEN=D2
  117.       IF (DEN) 90,80,90
  118.  80   CA=0.
  119.       CG=-A(J,K)/A(K,K)
  120.       GO TO 100
  121.  90   CA=AKK/DEN
  122.       CG=-AJJ/DEN
  123. C
  124. C     WE PERFORM THE GENERALIZED ROTATION
  125.  100  IF (N-2) 101,180,101
  126.  101  JP1=J+1
  127.       JM1=J-1
  128.       KP1=K+1
  129.       KM1=K-1
  130. C
  131.       IF (JM1-1) 120,110,110
  132.  110  DO 105 I=1,JM1
  133.       AJ=A(I,J)
  134.       BJ=B(I,J)
  135.       AK=A(I,K)
  136.       BK=B(I,K)
  137.       A(I,J)=AJ+CG*AK
  138.       B(I,J)=BJ+CG*BK
  139.       A(I,K)=AK+CA*AJ
  140.  105  B(I,K)=BK+CA*BJ
  141. C
  142.  120  IF (KP1-N) 130,130,140
  143.  130  DO 125 I=KP1,N
  144.       AJ=A(J,I)
  145.       BJ=B(J,I)
  146.       AK=A(K,I)
  147.       BK=B(K,I)
  148.       A(J,I)=AJ+CG*AK
  149.       B(J,I)=BJ+CG*BK
  150.       A(K,I)=AK+CA*AJ
  151.  125  B(K,I)=BK+CA*BJ
  152. C
  153.  140  IF (JP1-KM1) 150,150,180
  154.  150  DO 160 I=JP1,KM1
  155.       AJ=A(J,I)
  156.       BJ=B(J,I)
  157.       AK=A(I,K)
  158.       BK=B(I,K)
  159.       A(J,I)=AJ+CG*AK
  160.       B(J,I)=BJ+CG*BK
  161.       A(I,K)=AK+CA*AJ
  162.  160  B(I,K)=BK+CA*BJ
  163.  180  AK=A(K,K)
  164.       BK=B(K,K)
  165.       A(K,K)=AK+2*CA*A(J,K)+CA*CA*A(J,J)
  166.       B(K,K)=BK+2*CA*B(J,K)+CA*CA*B(J,J)
  167.       A(J,J)=A(J,J)+2*CG*A(J,K)+CG*CG*AK
  168.       B(J,J)=B(J,J)+2*CG*B(J,K)+CG*CG*BK
  169.       A(J,K)=0.0
  170.       B(J,K)=0.0
  171. C
  172. C     UPDATE EIGENVECTORS
  173.       DO 190 I=1,N
  174.       XJ=X(I,J)
  175.       XK=X(I,K)
  176.       X(I,J)=XJ+CG*XK
  177.  190  X(I,K)=XK+CA*XJ
  178. C
  179.  210  CONTINUE
  180. C
  181.       DO 220 I=1,N
  182.  220  EIGV(I)=A(I,I)/B(I,I)
  183.       IF(IFPR.EQ.0) GO TO 227
  184.       WRITE (IOUT,1005)
  185.       WRITE (IOUT,1002) (EIGV(I),I=1,N)
  186.  227  CONTINUE
  187. C
  188. C     CHECK FOR CONVERGENCE
  189.       DO 240 I=1,N
  190.       TOL=RTOL*D(I)
  191.       DIF=DABS(EIGV(I)-D(I))
  192.       IF(DIF.GT.TOL) GO TO 300
  193.  240  CONTINUE
  194. C
  195. C     CHECK IF ALL OFF-DIAG ELEMENTS ARE SATISFACTORILY SMALL
  196.       EPS=RTOL**2
  197.       DO 260 J=1,NR
  198.       JJ=J+1
  199.       DO 260 K=JJ,N
  200.       TT=A(J,K)*A(J,K)
  201.       TB=A(J,J)*A(K,K)
  202.       EPSA=DABS(TT/TB)
  203.       TT=B(J,K)*B(J,K)
  204.       TB=B(J,J)*B(K,K)
  205.       EPSB=TT/TB
  206.       IF ((EPSA.LT.EPS).AND.(EPSB.LT.EPS)) GO TO 260
  207.       GO TO 300
  208.  260  CONTINUE
  209. C
  210.       DO 310 I=1,N
  211.       DO 310 J=I,N
  212.       B(J,I)=B(I,J)
  213.  310  A(J,I)=A(I,J)
  214.       RETURN
  215. C
  216.  300  DO 320 I=1,N
  217.  320  D(I)=EIGV(I)
  218.       IF (NSWEEP.LT.NSMAX) GO TO 40
  219.       DO 330 I=1,N
  220.       DO 330 J=I,N
  221.       B(J,I)=B(I,J)
  222.  330  A(J,I)=A(I,J)
  223.       RETURN
  224. C
  225.  1000 FORMAT (27H0SWEEP NUMBER IN *JACOBI* =, I4)
  226.  1002 FORMAT (1H ,12E11.4)
  227.  1004 FORMAT (37H0***ERROR   SOLUTION STOP IN *JACOBI*, / 12X,
  228.      1       8HCHECK = , E20.14 / 1X)
  229.  1005 FORMAT (36H0CURRENT EIGENVALUES IN *JACOBI* ARE, / 1X)
  230.  2020 FORMAT ('0**ERROR** MATRIX NOT POSITIVE DEFINITE IN JACOBI'/)
  231. C
  232.       END
  233.